clear all
tempfile tempsave
set seed ${seed}

do "${scripts}_output.do"
do "${scripts}_texgraph.do"

/*Macros*/
local data 		"National Longitudinal Study of Adolescent to Adult Health"
local source 	"Own calculations."
local signif 	"Significance levels:" ///
				"* \(p<0.10\), ** \(p<0.05\), *** \(p<0.01\)."
local error 	"Standard errors are clustered at the school level."
local controls 	"\emph{Child Controls}:" ///
				"Firstborn dummy, linear birth cohort trend (in months) by gender," ///
				"20 principal components of the full matrix of genetic data." ///
				"\emph{Family Controls}:" ///
				"Age of mother at birth, years of education of both mother and father," ///
				"average potential wages of both mother and father," ///
				"the standard deviation of potential wages of both mother and father, dummies for non-US born mothers and fathers," ///
				"a dummy for Christian religion, state fixed effects."
local controlfct	"\emph{Control Function}: Share white, share single mothers, maternal education (average),"  ///
					"share females, share migrants. All control function variables are calculated"  ///
					"as leave-cohort-out school averages."
local standard  "All right-hand side variables are standardized on the estimation sample ($\mu=0$, $\sigma=1$)."
local notesize scriptsize

local format pdf
set scheme white_tableau
graph set window fontface "Palatino Linotype"
colorpalette viridis, nograph
local c1="`r(p1)'"
local c2="`r(p5)'"
local c3="`r(p10)'"
local c4="`r(p15)'"

// -----------------------------------------------------
// Figure S.1
// -----------------------------------------------------
local name figS1

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace
xtile sesbelsky_cat=sesbelsky, nquantiles(10)
keep if inrange(pgs_edu,-4,4)

joyplot pgs_edu, by(sesbelsky_cat) yrev overlap(10) alpha(50) bwidth(0.5) lc(white%100) lw(0.4) ///
yl ylpattern(dot) ylabsize(3) ///
xtitle("PGI{superscript:EA}") ytitle("Decile of Family SES")
graph export "${graphpath}`name'.`format'", replace

local title "Distribution of \PGS by Family SES"
local desc  "This figure shows unconditional \PGS distribution by deciles of" ///
			"family SES.  Density distributions are smoothed"  ///
			"using the Epanechnikov kernel function with a bandwidth of 0.5."
texgraph, 	title(`title') name(`name') replace ///
		width(1) file(`format') notesize(`notesize')  ///
		note_paper("`source' `desc'") data(`data')  ///
		inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure S.2
// -----------------------------------------------------
local name figS2

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

qui: sensemakr eduyears f_over1 ///
$controls i.w1state, ///
treat(f_over1) benchmark(edufather) tcontourplot contourplot clim(0 0.5) kd(2(2)10) suppress

matrix A=e(contourgrid)
svmat A
matrix B=e(tcontourgrid)
svmat B
matrix C=e(bounds)
svmat C
rename 	(A1 A2 A3 B3 C1 C3 C4 C5 C7) ///
		(c_treat c_outc c_coef c_t m m_treat m_outc m_coef m_t)
keep c_treat c_outc c_coef c_t m m_treat m_outc m_coef m_t

replace c_t=2*ttail(e(dof),abs(c_t))
replace m_t=2*ttail(e(dof),abs(m_t))

local maxdim 0.06
local maxax= `maxdim'+0.005
local cuts `cuts' 0.01 0.05 0.1 1
local original="`: di string(e(treat_coef),"%9.3f")'"

tostring m, replace
gen marker = m+"x"
replace marker=marker + " (Paternal Education)" if marker=="10x"

graph twoway ///
(contourline c_t c_outc c_treat if c_coef>0 & c_treat<`maxdim' & c_outc<`maxdim', ccuts(`cuts')) ///
(scatter m_outc m_treat, mlabel(marker) mcolor("`c2'") mlabcolor("`c2'")) ///
(scatteri 0 0, mcolor("`c1'") msym(diamond)), ///
xscale(range(0(0.02)`maxax')) xlabel(0(0.02)`maxax') ///
yscale(range(0(0.02)`maxax')) ylabel(0(0.02)`maxax') ///
xtitle("Partial R{superscript:2} of Unobserved Confounders with Q{superscript:S}") ytitle("Partial R{superscript:2} of Unobserved Confounders with Years of Educ.") ///
legend(off) ///
text(0.00 0.000 " `original' (Baseline)", place(e) size(small) color(`c1')) ///
text(0.006 0.006 "p=0.05", place(ne) size(small) color(gs8)) ///
text(0.012 0.012 "p=0.10", place(ne) size(small) color(gs8)) ///
text(0.039 0.039 "{&beta}=0", place(ne) size(small) color(gs8))

graph export "${graphpath}`name'.`format'", replace

local title "Sensitivity of \IQuali to Unobserved Confounders"
local desc  "This figure shows the sensitivity of the point estimate for the impact of" ///
			"\IQuali on years of education to unobserved confounding variables." ///
			"The diamond shows baseline estimates from Figure \ref{fig: fig1}." ///
			"\citet{Cinelli2020} show that omitted variable bias is a function of two partial $ R^2 $:" ///
			"(i) the partial $ R^2 $ of unobserved confounders with the outcome variable conditional on controls (y-axis)," ///
			"(ii) the partial $ R^2 $ of unobserved confounders with the treatment variable conditional on controls (x-axis)." ///
			"We plot three contour lines showing for which combinations of (i) and (ii), the coefficient of \IQuali" ///
			"would drop below a significance level of 5\%, 10\%, and when it would flip signs." ///
			"For comparison, we also plot adjusted treatment effects of \IQuali assuming that we would control for" ///
			"an unobserved confounder with values of (i) and (ii) that are a (2x-, 4x-, 6x-, 8x-, 10x-) multiple of paternal education (measured in years)." ///
			"This comparison suggests that an observed confounder even 10 times as strong as paternal education could not drive" ///
			"the coefficient of \IQuali below a significance level of 5\%." ///
			"This result is mainly driven by the low partial $ R^2 $ of paternal education with \IQuali after partialling out our controls."
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize')  ///
			note_paper("`source' `desc' `controls' `controlfct' `standard' `error'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure S.3
// -----------------------------------------------------
local name figS3

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

local quantiles 10

reghdfe pgs_edu  $controls [pw=w0], absorb(w1state) cluster(HsId) resid
predict res_pgs, residuals
reghdfe f_over1 $controls  [pw=w0], absorb(w1state) cluster(HsId) resid
predict res_f_over1, residuals


xtile f_over1_cat=res_f_over1, nquantiles(10)

local total=0
local ten=0
local five=0
local one=0
forvalues i=1/10{
	local k=`i'+1
	forvalues j=`k'/10{
		ksmirnov res_pgs if f_over1_cat==`i' | f_over1_cat==`j', by(f_over1_cat)

		local total=`total'+1

		if `r(p)'<=0.1 & `r(p)'>0.05{
			local ten=`ten'+1
		}
		if `r(p)'<=0.05 & `r(p)'>0.01{
			local five=`five'+1
		}
		if `r(p)'<=0.01{
			local one=`one'+1
		}
	}
}

di `total'
di `ten'
di `five'
di `one'

keep if inrange(pgs_edu,-4,4)

joyplot res_pgs, by(f_over1_cat) yrev overlap(10) alpha(50) bwidth(0.5) lc(white%100) lw(0.4) ///
yl ylpattern(dot) ylabsize(3) ///
xtitle("PGI{superscript:EA}") ytitle("Decile of Q{superscript:S}")
graph export "${graphpath}`name'.`format'", replace

local title "Distribution of \PGS by \IQuali"
local desc  "This figure shows \PGS distribution by deciles of" ///
			"\IQuali\negspace after residualizing \PGS and \IQuali by the full set of control variables." ///
			"Density distributions are smoothed"  ///
			"using the Epanechnikov kernel function with a bandwidth of 0.5."
texgraph, 	title(`title') name(`name') replace ///
		width(1) file(`format') notesize(`notesize')  ///
		note_paper("`source' `desc'") data(`data')  ///
		inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure S.4
// -----------------------------------------------------
local name figS4

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

foreach x in hs coll4y postgrad pvt pat risk{
	reghdfe `x' pred [aw=w0], absorb(w1state) cluster(HsId)
	local b=string(_b[pred],"%9.3f")
	local se=string(_se[pred],"%9.3f")
	local N=string(e(N),"%9.0fc")
	local R=string(e(r2),"%9.3f")

	if "`x'"=="hs"{
			local gtitle "High School Degree"
			local h 1.1
	}
	if "`x'"=="coll4y"{
			local gtitle "4-year College Degree"
			local h 1.0
	}
	if "`x'"=="postgrad"{
			local gtitle "Postgraduate Degree"
			local h 0.6
	}
	if "`x'"=="pvt"{
			local gtitle "Picture Vocabulary Test"
			local h 1.0
	}
	if "`x'"=="pat"{
			local gtitle "Patience"
			local h 0.5
	}
	if "`x'"=="risk"{
			local gtitle "Risk Aversion"
			local h 0.2
	}

	binscatter `x' pred[aw=w0], ///
	lcolor(gs0) mcolor(edkblue) msymbols(O) ///
	xtitle(" ") ytitle("") ///
	title("`gtitle'", size(medsmall) color(gs0)) ///
	text(`h' 12.0 "R{superscript:2}=`R'", place(e) color(gs0) just(left) margin(l+1 t+1 b+1 r+1) fcolor(gs0%0) size(small))
	graph save "${graphpath}pred_`x'.gph", replace
}
graph combine 	"${graphpath}pred_hs" "${graphpath}pred_coll4y" "${graphpath}pred_postgrad" ///
				"${graphpath}pred_pvt" "${graphpath}pred_pat" "${graphpath}pred_risk"  ///
				, cols(3) xcommon b1("Predicted Education (in Years)") imargin(0 0 0 0)
graph display , xsize(4)
graph export "${graphpath}`name'.`format'", replace

local title "Association of Educational Attainment and Skills with Predicted Education"
local desc  "This figure shows the association of various meausures for educational" ///
			"attainment and (non-)cognitive skills with predicted years of education." ///
			"We bin scatterplots using 20 quantiles of the variable of interest." ///
			"Black lines are fitted from linear regressions of the variable of interest on predicted education." ///
			"Predicted education is calculated from a regression of" ///
			"completed years of education on all \emph{Child Controls} and \emph{Family Controls}." ///
			"Outcomes are measured in Waves 3, 4 and 5."
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize')  ///
			note_paper("`source' `desc' `controls'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure S.5
// -----------------------------------------------------
local name figS5

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

cap postfile figs_5 hs be cu cl using "${temp}figs_5.dta", replace

local m=0
reghdfe eduyears c.pgs_edu##(c.f_over1) ///
$controls [pw=w0], absorb(w1state) cluster(HsId) keepmata
local b=_b[c.pgs_edu#c.f_over1]
local cu=_b[c.pgs_edu#c.f_over1]+_se[c.pgs_edu#c.f_over1]*invnorm(0.95)
local cl=_b[c.pgs_edu#c.f_over1]-_se[c.pgs_edu#c.f_over1]*invnorm(0.95)
post figs_5 (`m') (`b') (`cu') (`cl')
local ++m

glevelsof HsId, local(schools)
foreach x in `schools'{
	preserve
	drop if HsId==`x'
	reghdfe eduyears c.pgs_edu##(c.f_over1) ///
	$controls [pw=w0], absorb(w1state) cluster(HsId) keepmata
	local b=_b[c.pgs_edu#c.f_over1]
	local cu=_b[c.pgs_edu#c.f_over1]+_se[c.pgs_edu#c.f_over1]*invnorm(0.95)
	local cl=_b[c.pgs_edu#c.f_over1]-_se[c.pgs_edu#c.f_over1]*invnorm(0.95)
	post figs_5 (`m') (`b') (`cu') (`cl')
	restore
	local ++m
}
postclose figs_5

use "${temp}figs_5.dta", clear
lab def out 0 "Baseline"
lab val hs out

graph twoway ///
(rcap cu cl hs if hs==0, lcolor("`c1'"%30))  ///
(scatter be hs if hs==0, mcolor("`c1'") msymbol(O) msize(large))  ///
(rcap cu cl hs if hs!=0, lcolor("`c3'"%30) lpattern(shortdash_dot))  ///
(scatter be hs if hs!=0, mcolor("`c3'") msymbol(s) msize(small)) ///
,  ///
yline(0, lpattern(dash))  ///
ytitle("Coefficient of PGI{superscript:EA} x Q{superscript:S}")  ///
ylabel(-0.125(0.025)0, format(%9.3fc)) ///
xtitle("ID of Dropped High School") ///
xlabel(0(10)73, valuelabels) ///
legend(off)
graph export "${graphpath}`name'.`format'", replace

local title "Sensitivity to Outlier Schools"
local desc  "This figure shows point estimates and 90\% confidence bands of the" ///
			"interaction of \PGS with \IQuali\negspace, and its association" ///
			"with years of education. Each estimate is derived from a subsample of the data" ///
			"in which we drop one High School, respectively." ///
			"Estimates follow the specification of equation (\ref{eq: estim_educ})."
texgraph, 	title(`title') name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc' `controls' `controlfct' `standard' `error'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure S.6
// -----------------------------------------------------
local name figS6

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

cap postfile figs_6 bound be cu cl using "${temp}figs_6.dta", replace
gen eduyears_b=eduyears
forvalues i=20(-1)8{
	replace eduyears_b=`i' if eduyears_b>`i'
	reghdfe eduyears_b c.pgs_edu##(c.f_over1) ///
	$controls [pw=w0], absorb(w1state) cluster(HsId) keepmata

	local b=_b[c.pgs_edu#c.f_over1]
	local cu=_b[c.pgs_edu#c.f_over1]+_se[c.pgs_edu#c.f_over1]*invnorm(0.95)
	local cl=_b[c.pgs_edu#c.f_over1]-_se[c.pgs_edu#c.f_over1]*invnorm(0.95)
	post figs_6 (`i') (`b') (`cu') (`cl')

}
postclose figs_6

use "${temp}figs_6.dta", clear
lab def ceil 20 "Baseline"
lab val bound ceil

graph twoway ///
(rcap cu cl bound if bound==20, lcolor("`c1'"%30))  ///
(scatter be bound if bound==20, mcolor("`c1'") msymbol(O) msize(large))  ///
(rcap cu cl bound if bound!=20, lcolor("`c3'"%30) lpattern(shortdash_dot))  ///
(scatter be bound if bound!=20, mcolor("`c3'") msymbol(s) msize(small)) ///
,  ///
yline(0, lpattern(dash))  ///
ytitle("Coefficient of PGI{superscript:EA} x Q{superscript:S}")  ///
ylabel(-0.125(0.025)0, format(%9.3fc)) ///
xtitle("Censoring Level (Years of Education)") ///
xlabel(20(1)8, valuelabels) xscale(reverse) ///
legend(off)
graph export "${graphpath}`name'.`format'", replace

local title "Sensitivity to Ceiling Effects"
local desc  "This figure shows point estimates and 90\% confidence bands of the" ///
			"interaction of \PGS and \IQuali\negspace, and its association" ///
			"with years of education. Each estimate is derived from the full sample while" ///
			"censoring the outcome variable at different levels from above." ///
			"Estimates follow the specification of equation (\ref{eq: estim_educ})."
texgraph, 	title(`title') name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc' `controls' `controlfct' `standard' `error'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs}) 


// -----------------------------------------------------
// Figure S.7
// -----------------------------------------------------
local name figS7

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

reghdfe pvt (c.pgs_edu)##(c.f_over1) ///
$controls [pw=w0], absorb(w1state) cluster(HsId)
margins, at(pgs_edu=(-4.0(.5)3.0) f_over1 =(-2.5(.5)1.5)) saving("${temp}figs_7", replace)

use "${temp}figs_7", clear
colorpalette viridis, nograph
graph twoway ///
(contour _margin _at1 _at2, levels(10) ///
ccolors( "`r(p1)'" "`r(p2)'"  "`r(p3)'"  "`r(p4)'"  "`r(p6)'" ///
		 "`r(p8)'"  "`r(p10)'"  "`r(p12)'"  "`r(p14)'" "`r(p15)'")) ///
, ///
xtitle("Q{superscript:S}", size(large)) xsize(8) ///
ytitle("PGI{superscript:EA}", size(large)) yscale(range(-4(.5)3.0)) ylabel(-4(1)3.0) xscale(range(-2.5(1)1.5)) xlabel(-2.5(1)1.5) ///
ztitle("PVT" "(Prediction)") zlabel(, format(%9.1f)) ///
clegend(position(4) bmargin(zero))
graph export "${graphpath}`name'.`format'", replace

local title "Association of Picture Vocabulary Test with \PGS by \IQuali"
local desc  "This figure shows predictions of the PVT" ///
			"by \PGS and \IQuali cell." ///
			"Predictions are calculated using the model estimated in column (1) in Panel (a)" ///
			"of Table \ref{tab: tab7} while allowing for non-linear effects over the range of -4.0(0.5)3.0 of \PGS and -2.5(0.5)1.5 of \IQuali\negspace."
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})


// -----------------------------------------------------
// Figure S.8
// -----------------------------------------------------
local name figS8

local reps=10000
use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

permute f_over1 _b[c.pgs_edu#c.f_over1] _se[c.pgs_edu#c.f_over1], reps(`reps') rseed($seed) saving("${temp}permute.dta", replace): ///
reghdfe eduyears c.pgs_edu##(c.f_over1)  ///
$controls [pw=w0], absorb(w1state) cluster(HsId)

reghdfe eduyears c.pgs_edu##(c.f_over1)  ///
$controls [pw=w0], absorb(w1state) cluster(HsId)
local tbase=_b[pgs_edu#f_over1]/_se[pgs_edu#f_over1]

use "${temp}permute.dta", replace
gen t=_pm_1/_pm_2
local var t
hashsort `var'
gen order=_n
sum `var' if order==(1/100*`reps')
local i1=r(mean)
sum `var' if order==(2.5/100*`reps')
local i5=r(mean)
sum `var' if order==(5/100*`reps')
local i10=r(mean)

graph twoway (hist `var', fcolor(white%0) lcolor(tile) percent ///
			xline(`tbase', lcolor(maroon) lwidth(thick)) ///
			xline(`i1', lcolor(blue%80) lpattern(solid)) ///
			xline(`i5', lcolor(blue%50) lpattern(solid)) ///
			xline(`i10', lcolor(blue%30) lpattern(solid))), ///
			xtitle("t-value") ///
			text( 4 -2.40 "1%", color(blue%80) orientation(vertical)) ///
			text( 4 -1.90 "5%", color(blue%50) orientation(vertical)) ///
			text( 4 -1.50 "10%", color(blue%30) orientation(vertical)) ///
			text( 4 -2.80 "Baseline estimate", color(maroon) orientation(vertical))


graph export "${graphpath}`name'.`format'", replace

local title "Permutation test for placebo assignments of \IQuali"
local desc  "This figure shows the frequency distribution of $ t $-statistics for the" ///
			"interaction of \PGS and \IQuali under 10,000 permutations of \IQuali\negspace."
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc' `controls' `controlfct' `standard' `error'")  data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure S.9
// -----------------------------------------------------
local name figS9


use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace


foreach x in health_sub health_obj pvt risk pat intel consc extra agree neuro{
	reghdfe eduyears `x' $controls [pw=w0], absorb(w1state) cluster(HsId) keepmata
	local N=string(e(N),"%9.0fc")
	local R=string(e(r2),"%9.3f")
	local b=string(_b[`x'],"%9.3f")
	local se=string(_se[`x'],"%9.3f")

	if "`x'"=="pvt"{
			local gtitle "Picture Vocabulary Test"
	}
	if "`x'"=="risk"{
			local gtitle "Risk Aversion"
	}
	if "`x'"=="pat"{
			local gtitle "Patience"
	}
	if "`x'"=="health_obj"{
			local gtitle "Health (Obj.)"
	}
	if "`x'"=="health_sub"{
			local gtitle "Health (Subj.)"
	}
	if "`x'"=="intel"{
			local gtitle "Openness"
	}
	if "`x'"=="consc"{
			local gtitle "Conscientiousness"
	}
	if "`x'"=="extra"{
			local gtitle "Extraversion"
	}
	if "`x'"=="agree"{
			local gtitle "Agreeableness"
	}
	if "`x'"=="neuro"{
			local gtitle "Neuroticism"
	}

	binscatter eduyears `x', controls($controls) absorb(w1state) ///
	lcolor(gs0) mcolor(edkblue) msymbols(O) ///
	xtitle(" ") ytitle("") ylabel( ,format(%9.1fc)) ///
	title("`gtitle'", size(medsmall) color(gs0)) ///
	text(15.9 -3.1  "Slope: `b'*** (`se')", place(e) color(gs0) just(left) margin(l+1 t+1 b+1 r+1) fcolor(gs0%0) size(small))
	graph save "${graphpath}mech_`x'.gph", replace
}

graph combine 	"${graphpath}mech_pvt" "${graphpath}mech_risk" "${graphpath}mech_pat" ///
				"${graphpath}mech_health_sub" "${graphpath}mech_intel" "${graphpath}mech_neuro"  ///
				, cols(3) ycommon xcommon l1("Education (in Years)") imargin(0 0 0 0)
graph display , xsize(4)
graph export "${graphpath}`name'.`format'", replace

local title "Association of Educational Attainment with (Non-)Cognitive Skills"
local desc  "This figure shows the association of educational" ///
			"attainment and various measures of (non-)cognitive skills." ///
			"We bin scatterplots using 20 quantiles of the variable of interest." ///
			"Black lines are fitted from linear regressions of educational attainment on" ///
			"the variable of interest and the full set of controls." 
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize')  ///
			note_paper("`source' `desc' `controls' `controlfct'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure S.10/S.11/S.12
// -----------------------------------------------------

// Simulation Bias
clear
set obs 4034
local iterations=50
gen x1=rnormal(0,1)
gen x2=rnormal(0,1)
gen z=.
gen x1x2=.
gen x1z=.
gen x2z=.
gen y=.

//runtime of simulation approx. 11h
cap postfile bias w1 w2 d1 d2 d3 cov1 cov2 iteration b1 b2 b3 r2 using "${temp}figs_10.dta", replace

qui forvalues w1=0.0(0.05)0.55{
	forvalues w2=0.0(0.05)0.55{

		noi di "Cov1=`w1' Cov2=`w2'"

		replace z=`w1'*x1+`w2'*x2+sqrt(1-`w1'^2-`w2'^2)*rnormal(0,1)
		replace x1x2=x1*x2
		replace x1z=x1*z
		replace x2z=x2*z

		// local vars x1 x2 z x1x2 x1z x2z
		// sum `vars'
		// corr `vars', cov
		corr x1 z, cov
		local cov1=r(cov_12)
		corr x2 z, cov
		local cov2=r(cov_12)

		forvalues d1=0.2(0.1)0.2{
			forvalues d2=0.0(0.01)0.21{
				forvalues d3=-0.21(0.01)0.1{
					forvalues i=1/`iterations'{
						replace y=0.4*x1+0.2*x2-0.1*x1x2+`d1'*z+`d2'*x1z+`d3'*x2z+rnormal(0,1)

						reg y x1 x2 x1x2

						post bias (`w1') (`w2') (`d1') (`d2') (`d3') (`cov1') (`cov2') (`i') (_b[x1]) (_b[x2]) (_b[x1x2]) (`e(r2)')
					}
				}
			}
		}
	}
}
cap postclose bias

local name figS10

use using "${temp}figs_10.dta", clear
gcollapse b1 b2 b3 r2 cov1 cov2, by(w1 w2 d1 d2 d3)

gen prod1=d2*w2
gen prod2=d3*w1
gen test=(prod1+prod2>0)
sum test

twoway (hist b3 if test==1, percent lcolor(gs12) fcolor(gs12%20)) ///
	   (hist b3 if test==0, percent fcolor(white%0) lcolor(tile) xline(-0.1, lcolor(maroon) lwidth(thick))), ///
	   xtitle("Estimate of {&beta}{subscript:3}", size(normal)) xlabel(-.3(0.1)0.1) ///
	   legend(off) text(8 -0.075 "{&beta}{subscript:3}=-0.10", color(maroon) orientation(horizontal))
graph export "${graphpath}`name'.`format'", replace

local title "Simulation for $\beta_3=-0.10$"
local desc  "This figure presents estimates of $\beta_3$ under different assumptions about $\delta_2, \delta_3, \sigma^2_{x_1z}, \sigma^2_{x_2z}$." ///
			"The simulation is based on the following data generating process:" ///
			"$ y_i = 0.4 x_{1i} + 0.2 x_{2i} -0.1 (x_{1i} \cdot x_{2i}) + 0.2 z_i+ \delta_2 (x_{1i}\cdot z_i) + \delta_3 (x_{2i} \cdot z_i)  + \epsilon_i $," ///
			"where $\sigma^2_{x_1x_2}=0$, $\sigma^2_{x_1z}=\sigma^2_{x_2z}\in [0.0(0.05)0.50]$," ///
			"$\delta_2 \in [0.0(0.01)0.20]$, $\delta_3 \in [-0.20(0.01)0.00]$." ///
			"We run 50 iterations for each combination of parameters." ///
			"Gray-shaded bars represent estimates for parameter combinations such that $|(i)|<|(ii)|$." ///
			"Dark-lined bars represent estimates for parameter combinations such that $|(i)|>|(ii)|$." 
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc'")  data("Simulated data")  ///
			inpath(${graphpath}) path_paper(${paper_graphs})


local name figS11
use using "${temp}figs_10.dta", clear
gcollapse b1 b2 b3 r2 cov1 cov2, by(w1 w2 d1 d2 d3)

replace w1=round(w1*100)

forvalues i=5(5)50{

	sum b1 if w1==`i'
	local gennurt=floor(100*((r(mean)-0.4)/0.4))
	di `gennurt'

	hist b3 if w1==`i', percent fcolor(white%0) lcolor(tile) xline(-0.1, lcolor(maroon) lwidth(thick)) ///
	title("Genetic nurture: `gennurt' %", size(normal)) xlabel(-.3(0.1)0.1) ///
	xtitle("", size(normal))
	graph save "${graphpath}bias_2_`i'", replace
}
graph combine 	"${graphpath}bias_2_50" "${graphpath}bias_2_45" "${graphpath}bias_2_40" "${graphpath}bias_2_35" "${graphpath}bias_2_30" ///
				"${graphpath}bias_2_25" "${graphpath}bias_2_20" "${graphpath}bias_2_15" "${graphpath}bias_2_10" "${graphpath}bias_2_5" ///
				, cols(5) ycommon xcommon b1("Estimate of {&beta}{subscript:3}") imargin(0 0 0 0)
graph display , xsize(8)
graph export "${graphpath}`name'.`format'", replace

local title "Simulation for $\beta_3=-0.10$ under decreasing $\sigma^2_{x_1z}$"
local desc  "This figure presents estimates of $\beta_3$ under different assumptions about $\delta_2, \delta_3, \sigma^2_{x_1z}, \sigma^2_{x_2z}$." ///
			"The simulation is based on the following data generating process:" ///
			"$ y_i = 0.4 x_{1i} + 0.2 x_{2i} -0.1 (x_{1i} \cdot x_{2i}) + 0.2 z_i+ \delta_2 (x_{1i}\cdot z_i) + \delta_3 (x_{2i} \cdot z_i)  + \epsilon_i $," ///
			"where $\sigma^2_{x_1x_2}=0$, $\sigma^2_{x_1z}=\sigma^2_{x_2z}\in [0.0(0.05)0.50]$," ///
			"$\delta_2 \in [0.0(0.01)0.20]$, $\delta_3 \in [-0.20(0.01)0.00]$." ///
			"We run 50 iterations for each combination of parameters." ///
			"Each panel represent simulation results for a different specification of $\sigma^2_{x_1z}$" ///
			"over the interval $\sigma^2_{x_1z}=[0.05(0.05)0.50]$." ///
			"We denote the implied extent of confounding through genetic nurture in the table header."
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc'") data("Simulated data")  ///
			inpath(${graphpath}) path_paper(${paper_graphs}) 


local name figS12
use if w1==0 & w2<=0.50 using "${temp}figs_10.dta", clear
gcollapse b1 b2 b3 r2, by(w1 w2 d1 d2 d3)
colorpalette viridis, nograph
graph twoway ///
(contour b3 d2 w2, ///
level(10)  ///
ccolors( "`r(p15)'" "`r(p14)'"  "`r(p13)'"  "`r(p12)'" "`r(p11)'" "`r(p9)'" ///
		 "`r(p7)'" "`r(p5)'" "`r(p3)'"  "`r(p1)'")) ///
, ///
xtitle("{&sigma}{superscript:2}{subscript:x{subscript:2}z}", size(large)) xsize(8) ///
ytitle("{&delta}{subscript:2}", size(large)) ylabel(, format(%9.2f)) ///
ztitle("Estimate of {&beta}{subscript:3}", size(large)) zlabel(, format(%9.2f)) ///
clegend(position(4) bmargin(zero))

graph export "${graphpath}`name'.`format'", replace

local title "Simulation for $\beta_3=-0.10$ with $\sigma^2_{x_1z}=0$"
local desc  "This figure presents estimates of $\beta_3$ under different assumptions about $\delta_2, \delta_3, \sigma^2_{x_2z}$," ///
			"while fixing $\sigma^2_{x_1z}=0$." ///
			"The simulation is based on the following data generating process:" ///
			"$ y_i = 0.4 x_{1i} + 0.2 x_{2i} -0.1 (x_{1i} \cdot x_{2i}) + 0.2 z_i+ \delta_2 (x_{1i}\cdot z_i) + \delta_3 (x_{2i} \cdot z_i)  + \epsilon_i $," ///
			"where $\sigma^2_{x_1x_2}=0$, $\sigma^2_{x_1z}=0, \sigma^2_{x_2z}\in [0.0(0.05)0.50]$," ///
			"$\delta_2 \in [0.0(0.01)0.20]$, $\delta_3 \in [-0.20(0.01)0.00]$." ///
			"We run 50 iterations for each combination of parameters."
texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc'") data("Simulated data") ///
			inpath(${graphpath}) path_paper(${paper_graphs}) 


*delete temporary files
local list : dir "${graphpath}"  files "*.gph"
di `list'
foreach f of local list {
    erase "${graphpath}`f'"
}
